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Abstract 

Decision making at a cellular level determines different fates for isogenic cells. However, it is not yet clear how rational 
decisions are encoded in the genome, how they are transmitted to their offspring, and whether they evolve and become 
optimized throughout generations. In this paper, we use a game theoretic approach to explain how rational decisions are 
made in the presence of cooperators and competitors. Our results suggest the existence of an internal switch that operates 
as a biased coin. The biased coin is, in fact, a biochemical bistable network of interacting genes that can flip to one of its 
stable states in response to different environmental stimuli. We present a framework to describe how the positions of 
attractors in such a gene regulatory network correspond to the behavior of a rational player in a competing environment. 
We evaluate our model by considering lysis/lysogeny decision making of bacteriophage lambda in £ coli. 
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Introduction 

How do organisms adapt in order to survive in ever-changing 
environments? In this process, one important factor is mutation 
that affects the survival of organisms through proper genotypic 
changes. In a classical view, natural selection chooses the best-fit 
phenotype in the environment, and the next generation inherits 
the genotype associated with it. In a dynamic environment, 
genotype correction can be thought of an organism's strategy to 
survive. The fact that evolution occurs in an environment 
containing competing organisms, each playing its own optimized 
strategies, further complicates the picture. 

In terms of time, biological evolution is a process that stabilizes 
after several trial and errors. However, some environmental 
fluctuations happen so quickly that there is not enough time for the 
organism to adapt itself through appropriate mutations. Organ- 
isms may face challenges caused by fluctuations in extracellular 
conditions. Early studies focused on the relationship between these 
environmental fluctuations and genetic diversities [1-4]. Other 
studies have shown that phenotypic variation exists, not as a 
consequence of underlying heritable genetic variation, but as an 
independent way of adapting to an ever-changing environment 
[5-9] . They discovered that cells with the same genotype can play 
different strategies and exhibit different phenotypes even when 
they are living in an identical environment. 

From an extracellular point of view, phenotype variability as a 
risk-reducing strategy has been modeled by evolutionary game 



theory [10-16]. In this model, individuals with the same genotype 
compete for a longer life and more descendants by using various 
strategies. These strategies are interpreted as being different 
phenotypes. The fitness of an individual depends on the benefits 
and costs accrued by that individual in the presence of others. The 
pattern of phenotypic variations, which cannot be invaded by any 
alternative phenotypes, is described and predicted using the 
evolutionarily stable strategy (ESS) [10]. In ESS, the ratio of 
different phenotypes is considered as a mixed strategy which does 
not include changes at the level of the genotype. This concept has 
been successfully applied to explain why bacterial RNA-phage has 
different frequencies of Phi6 and PhiH2 phenotypes [11]. It has 
been shown that the fitness of phenotypes in RNA-phage 
generates a payoff matrix which is similar to the payoff matrix 
of the prisoner's dilemma problem. Furthermore, evolutionary 
game theory has provides an appropriate framework to learn 
important evolutionary phenomena such as altruistic behavior 
[12,13], the evolution of sex ratio [14], pathogen-host interaction 
[15], and the rate and quantity of ATP production in different 
pathways of ATP synthesis in yeast [16]. 

From an intracellular point of view, there exist several 
biochemical reactions underlying phenotypic variation. The study 
of these intracellular functions is associated with quantitative 
genetic analysis such as gene regulatory networks [17,18]. Gene 
regulatory networks play an important role in controlling the 
cellular behavior in varying environments. The structure and 
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features of gene regulatory networks are evolved as a result of 
adaptation to fluctuating environment [19]. 

A regulatory network of genes with positive and negative 
feedback loops creates a potential landscape with different 
attractor states and bifurcation points [20]. Moreover, it explains 
reshaping of the landscape based on alterations in network 
parameters. In the 1940s, Waddington provided a basis to explain 
how the cells of an organism evolve differently during embryonic 
development [21]. He introduced the term epigenetic landscape 
and portrayed it as a marble rolling down a mountain with 
different valleys that eventually comes to rest at the lowest point 
which represents the ultimate fate of the cells. The valleys of the 
landscape represent stable attractor states while the other less 
stable states represent transient states of the early embryonic or 
progenitor cells [20] . The location and the shape of the attractor 
states define the natural probabilities of different biological 
decisions. It is not entirely clear how a qualitative picture of a 
landscape can be quantified and how the structure of the 
landscape is encoded in the genome. 

We present a framework that combines the cooperative and 
competitive decision making of a living organism with its 
underlying intracellular gene regulatory network (see Figure 1). 
In this framework, game theoretic methods are applied to model 
the strategies of various living organisms. We show that the natural 
probabilities of organisms' decisions are fine-tuned to increase 
their chance of survival. Then, we argue that the location and the 
shape of the attractors in the Waddington landscape define the 
natural probabilities of different biological organisms, while the 
location and shape of attractors are characterized by the structure 
of the gene regulatory network. Altogether, we propose a 
framework that describes quantitatively how a gene regulatory 
network directs a cell to behave in a manner that is similar to that 
of a rational player in a game. This implies that the probability 
distribution of a rational decision, if we model a living organism as 
a cooperative and competitive decision maker, matches the 



probability distribution over stable states of its underlying gene 
regulatory network. 

Our framework is supported by experimental data from one of 
the well-studied biological cases of phenotypic variation, the 
infection of E. coli with bacteriophage lambda. After E. coli is 
infected by bacteriophage lambda, the virus chooses between 
lysogenic or lytic pathways [22-24]. In the lysogenic mode, the 
virus's genome is inserted into the bacterial genome and is 
replicated along with the bacterial genome; viral particles are not 
produced. In the lytic pathway, the viral genome is replicated 
independendy of the bacterial genome, viral particles are 
produced within the bacterial cell, the membrane of the host cell 
is lysed and the particles are released into the environment. Hence 
they can resume the infection cycle in other bacterial cells. 

The decision between lysogenic and lytic pathways contains a 
trade-off between safe maintenance of the viral genome within the 
bacterial host genome and increased bacteriophages production. 
Similar counter-intuitive situations exist where individuals sacrifice 
themselves and cooperate to win the overall game [25]. At the 
molecular level, the bacteriophage lambda gene regulatory 
network acts like a switch with two attractor states. This bistable 
network is the same as a biased coin that is falling on either side 
with different probabilities according to the environmental 
conditions. We demonstrate how the gene regulatory network 
directs phages to behave in a manner very comparable to what we 
expect from a rational player in a game. 

Description of the Model 

Representing living organisms as cooperative and competitive 
agents, game theory provides an elegant mathematical model to 
describe the behavior of a rational player. Consider an individual 
that makes a decision by selecting an action from a set of all 
possible actions A. The strategy of the individual can be seen as a 
probability distribution over all possible actions of the finite set A. 
In other words, the strategy of the individual can be defined by a 




Figure 1. The proposed framework. The framework that links the game theoretic perspective of decision making in a living organism and 

Waddington's perspective of the underlying gene sequence. 

doi:10.1371/journal.pone.0103569.g001 
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vector s : A-*[0, 1], where s(a) shows the probability of action a 
and ~^2 a e A s(a) = 1 . The strategy is pure if s(a) has only values 0 
or 1. It is mixed if s(a) accepts all values in the interval [0, 1]. The 
utility of an individual represents how efficient its behavior is with 
respect to the environment. In general, each individual's utility 
depends on the state of the environment and strategies of all other 
individuals in the game. Therefore, the utility of individual j is 
usually represented as a function iij : A A xf->R where JV is the 
set of size n of all individuals and £ is the set of all possible 
environmental states. In this paper, we consider a situation where 
the effects of other strategies are implicit, and assume the utility of 
an individual only depends on the environmental state, i.e., Uj{e) 
represents the utility of individual j when the state of the 
environment is e e £ . 

As a consequence, the utility of an action a e A depends on the 
state of the environment. However, an individual only detects 
internal signals through its biochemical network and thus predicts 
the environmental state. Let T be the set of all possible internal 
states and iel be an instance of it. Assume E and / are two 
random variables defining the environmental state and internal 
state respectively. An individual detects internal state i and infers 
the environmental state by determining posterior probabilities 
f(E = e\I = i), for each e e £ . To summarize the system dynamics, 
we express it in the following iterative steps. First, the environ- 
mental state e affects each individual's internal state. Second, each 
individual infers the environmental state based on its internal 
signals. Third, each individual estimates the expected utility of 
each action a based on its internal state. Fourth, each individual 
chooses its mixed strategy based on the expected utility of each 
action. At last, all individuals' strategies determine the future 
environmental state. It is worth mentioning that while individuals' 
strategies affect the environmental state, we do not aim to model 
these effects in this paper. 

From an intracellular point of view, the gene regulatory network 
is responsible for determining an organism's behavior. We 
represent the gene regulatory network Q as a set of m vertices 
V = {v\, V2, ■ ■ ■ , v m } and a set of weights W where Wy is the 
effect of gene v,- on gene v,-. A positive value for Wu means la- 
minates the expression of Vj and a negative value for Wjj means V; 
inhibits the expression of Vj. Every gene regulatory network 
characterizes a dynamical system which can be represented by a 
set of differential equations [26,27]. Let Vg be the set of all 
attractors of the dynamical system of gene regulatory network Q. A 
dynamical system in good conditions ends up in one of its attractor 
states as time goes to infinity and each attractor state induces a 
specific action. In this model, the action of a cell deterministically 
depends on its initial conditions. However, finding the exact initial 
conditions of a cell in a noisy environment is impractical. 
Therefore, we assign different probabilities to different attractor 
states by considering stochastic initial states. In particular, a 
probability distribution q : X>g— >[0, 1] is assigned to a dynamical 
system to describe the state of the system as times goes by. For 
every d e Vg the value q(d) shows the probability that a system 
with a random starting point ends up in attractor state d. 

Results 

Let us consider an organism which receives internal signals and 
selects between all its possible actions. The probability distribution 
s will be the strategy of a rational player in this situation. We 
demonstrate that for every internal, signal the weights of edges in 
the gene regulatory network are changed such that the probability 
distribution q of the corresponding dynamical system is almost the 
same as the mixed strategy s of the corresponding game. 



Gene regulatory network of the lysis/lysogeny decision 

After the E. coli has been infected with bacteriophage lambda, a 
decision is made between either the death of the host (lysis) or viral 
dormancy (lysogeny) [28,29] . Multiple internal and environmental 
signals arising during the decision making process affect the 
biochemical network and consequently, the fate of phages. These 
signals include the number of infecting phages, the metabolic state 
of the bacterium, the position of the infecting phages on the 
bacterium surface, and the bacterium size [22-24,30-32]. To 
model the decision making process of bacteriophage lambda, some 
parameters are more prominent in the final state of bacteriophage. 
The number of phages infecting a bacterium (multiplicity of 
infection; MOI) has long been known to affect the final decision of 
the bacteriophage [22,30-32]. In addition, recent results represent 
the volume of the infected bacterium as an important parameter in 
the decision-making process [23,24]. Recent results also show the 
relation between the concentration of phages in the host bacterium 
and the phage's final fate [23,24]. 

From intracellular perspective, the decision making of bacte- 
riophage lambda is linked to the competitive gene expression of cl 
and Cro (see Figure 2). The regulatory circuit chooses between 
two outcomes including cl and Cro genes. The expression of the 
cl gene will lead to a lysogeny decision, and Cro expression is 
followed by lysis [30-34]. These two genes are about 100 
nucleotides far apart in the genome, and the area between them 
consists of the promoter of Cro (Pr), the promoter of cl (Prm), 
and three operator sites ORi, OR2, and OR3. The operator sites 
ORi and OR2 are located in Pr, and operator site OR3 is located 
in Prm [35] (see Figure 2). Dimer CI2 has a strong binding affinity 
to ORi. Dimer Cr02 acts contrarily with the highest binding 
affinity to OR3. Binding of CI2 to ORi and OR2 inhibits Pr (down 
regulating Cro) and has a positive effect on Prm (auto-activating), 
while binding of Cr02 to the operator sites has a negative effect on 
the expression of the two genes. 

Rational decision making perspective 

After the infection of E. coli by bacteriophage lambda, a 
decision is made between lysis and lysogeny. A rational decision 
infers the environmental state (the average MOI), based on 
internal signals and increases the probability of lysogeny in a 
higher than expected average MOI. The host bacterium MOI, 
size, and concentration as internal signals define the phage fate. 
We calculate the probability distribution of estimated average 
MOI in different situations and describe the effect of host 
bacterium MOI, size and the concentration on it. Since the 
probability of each action, s(a), highly depends on estimated 
average MOI, a higher estimated average MOI results in a higher 
probability of lysogeny. 

Let f{E = X\I L = j£) be the probability that the average MOI 
equals X, given host bacterium size L and host bacterium MOI /.(. 
In other words, we are infering the enviromental state (the avarege 
MOI), with respect to the internal signals (host bacteriom size and 
MOI). The probability/^ = X\Il = £i) is calculated in Figure 3 for 
a fixed n and different host bacterium sizes. The estimated average 
MOI is highly correlated with the number of phages in a fixed size 
host bacterium. 

The probability f(E = X\lL = Lfi) is shown in Figure 3 for a 
specific host bacterium concentration and different host bacterium 
sizes. Figure 3 demonstrates the effect of the host bacterium size 
on the estimation of X for the same host bacterium concentration. 
The expected value of estimated average MOI remains un- 
changed for the same host bacterium concentrations. However, it 
has been shown that for the same concentration, the variance 
decreases as the size of the host bacterium increases. This implies 
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Lysogenic 



Figure 2. Regulatory network of bacteriophage lambda. Left: Genes and operators which are involved in lysogenic and lytic process of 

bacteriophage Lambda. Right: Simplified regulatory network of cl and Cro in bacteriophage lambda. 

doi:10.1371/journal.pone.0103569.g002 



that increasing the size of the host bacterium decreases the 
probability of having a high environmental MOI. In conclusion, 
we expect different reactions from a rational player in situations 
with the same host bacterium concentration but different host 
bacterium sizes. 

In Figure 4, we use the experimental data of [30-32] and [24] 
to analyze the phages' behavior. The experimental data of [30-32] 
measures the rate of lysogeny versus the average MOI. The same 
data has been used in [36]. The second experimental data set 
comes from [24]. This work measures both the probability of 
lysogeny based on the host bacterium MOI and also the rate of 
lysogeny based on the average MOI [24]. We calculate the 
expected probability s(a) for both the lysis and the lysogeny 
actions. For a given average MOI, we integrate over all possible 
host bacterium sizes and different MOIs to calculate the expected 
probability of lysogeny. The probability distribution of the host 
bacterium MOI is modeled by a Poisson distribution for a given 
average MOI. Figure 4A shows the expected probability of 
entering lysogeny as a function of the average MOI. This figure 
clearly illustrates that the behavior of phages is remarkably close to 
the behavior of rational players. We use the root-mean-square 
error (RMSE) to measure the accuracy of the proposed model and 
obtain an RMSE of 0.105. 

We use the experimental data of [24] to verify our results. The 
data in [24] estimates the probability of entering lysogeny based on 
the host bacterium MOI. We compute the expected probability 
s(a) for both actions lysis and lysogeny. For a given host bacterium 
MOI, we integrate over all sizes of the host bacterium to calculate 
the expected probability of lysogeny (see Figure 4). An RMSE of 



0.048 confirms that the behavior of a rational player is almost the 
same as the bacteriophages' behavior. 

Figure 4i? also shows the probability of lysogeny for different 
host bacterium MOIs. Phages tend to lysogenize as the MOI 
increases, since the estimation of average MOI is an increasing 
function of the host bacterium MOI. Moreover, Figure 4C 
indicates that the probability of lysogeny decreases as the size of 
the host bacterium increases. In fact, the growth of the host 
bacterium size reduces both the concentration of phages and the 
estimated average MOI. 

The probability of lysogeny in the different concentration of 
phages inside the host bacterium is illustrated in Figure 4Z). As we 
have discussed, the size of the host bacterium is one of the most 
important variables in the lysis/lysogeny decision [23]. Neverthe- 
less, in the same concentration of the host bacterium, the phages in 
the smaller host bacterium show more tendency to lysogenize [24] . 
Our model's assessment of phage behavior with an RMSE of 
0.077 matches the results of [24] as shown in Figure 4. Let us 
recall that the variance of average MOI estimation is higher for a 
small bacterium than a large bacterium (see Figure 3). Hence, in a 
fixed concentration, it is more likely to have a higher estimated 
avarage MOI in a smaller host bacterium. This means in the same 
concentration, the probability of lysogeny increases as the size of 
the host bacterium decreases. 

Intracellular perspective 

Analytical and numerical methods for solving differential 
equations will enable us to derive a set of steady state equations 
for a given system. To analyze this dynamic system, however, we 
would need more advanced methods. For this reason, we 




Figure 3. The probability distribution of estimated average MOI. (A) The host bacteria have the same MOI, but different sizes. (B) The host 

bacteria have the same phage concentration, but different sizes. 

doi:10.1371/journal.pone.0103569.g003 
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Figure 4. The probability of lysogeny as a function of the avarege phage input, the host bacterium MOI, and the host bacterium 
concentration. (A) The average phage input is the total number of phages divided by the total number of bacteria and it is shown in logarithmic 
scale. Red circles show phages behavior based on the experimental data of [30-32]. Green circles show phages behavior based on the experimental 
data of [24]. Blue curve shows our estimation of rational players reactions with a RMSE of 0.105. (B) The probability of entering lysogeny as a function 
of the host bacterium MOI. The host bacterium MOI is shown in logarithmic scale. Red circles show phages behavior base on the experiments of 
[24]. Blue curve shows the estimation of rational players reactions with an RMSE of 0.048. (C) The effect of the size of host bacterium on the 
probability of lysogeny. The results are displayed for the host bacterium sizes 0.7, 0.9, and 1.1. (D) The probability of lysogeny versus host bacterium 
concentration with different host bacterium MOI. Blue, red, and black lines represent the estimations based on our model for MOI = l, 2, and 3 
respectively. Blue, red, and black points show the experimental results of [24] for MOI = 1, 2, and 3 respectively. Blue plus shows the experimental 
result of [30-32] for MOI= 1. Note that our estimations have an RMSE of 0.077. 
doi:10.1371/journal.pone.0103569.g004 



approach the problem from a Waddington perspective. As shown 
in Figure 5, the Waddington landscape of the phages' decision- 
making varies with the number of phages infecting the same 
bacterium. In the case of a single phage infection, there is one 
attractor state with high Cro concentration that results in a lysis 
decision. Thus, the majority of the bacteria undergo lysis after a 
single infection (although there is a small probability for lysogeny 
due to the internal noise). A new attractor state with high cl 
concentration is distinguishable at MOI = 3. This results in a 
lysogeny decision. At MOI = 6, although there is still a vestige of 
the lysis attractor, but the deeper valley of the other attractor 
shows a higher likelihood of attracting the phages toward lysogeny. 
Video S 1 describing changes of attractor points in the Waddinton 
landscape, for the host bacterium MOI from 1 to 7, is provided in 
the supplementary material. 

To visualize the proposed landscape similar to the well-known 
drawing of Waddington envisioned in his book [21], the results of 
different MOI values are integrated in a 3D contour (Figure 6). 
According to the Waddington's original drawing, the vertical axis 
of the diagram represents the potential function and the back-to- 
front axis represents time [20]. For the lysis/lysogeny decision, we 
assign different values of MOI to the back-to-front axis. The left- 
to-right axis of Waddington's original diagram is the output of the 
system. It can be either the gene expression levels or phenotypic 
markers, whichever suitably distinguish between different cell 
types. The concentration level of the cl gene is chosen for this 
purpose in Figure 6. The attractors with low cl concentration 
represent the lysis state, and the attractor with high cl 
concentration represents the lysogeny state. 

The Waddington landscape of a host bacterium with different 
MOI is shown in Figure 6. In a single phage infection, (MOI = 1), 
the bacterium goes for lysis in the backside of the landscape. There 



is only one lysis valley with cl concentration close to 0. Due to the 
stochasticity of the system, there is a low probability for bacterium 
to undergo lysogeny after a single infection. As MOI increases, the 
single valley bifurcates and a new attractor for lysogeny appears in 
high cl concentrations. The front-side of the landscape shows what 
happens at MOI = 7. At this point, the channel for lysogeny is 
deeper (more stable) and wider (more probable) than the lysis. The 
high-potential barrier separating the two valleys represents the 
unstable states which the cells rarely pass. It justifies why the 
bacterium cannot easily change the outcome from lysis to lysogeny 
or vice versa. 

For every value of MOI from 0 to 7 with 0.1 steps, we 
performed 100000 simulation iterations to compute the probabil- 
ity of lysis or lysogeny decision-making (Figure 7). In each 
iteration, we pick a random starting point in the landscape for 
the initial expression values of cl and Cro, and track the trajectory 
of the cell towards different attractor states for a limited time, 
based on equation set 12. The commitment to either lysis or 
lysogeny is decided based on the final relative abundance of cl and 
Cro proteins (> 1 for lysogeny, < 1 for lysis). For every value of 
MOI, the lysis or lysogeny probabilities are defined as the ratio of 
the trajectories that ended in each corresponding state (Figure 7). 
We selected 100000 iterations since the probabilities do not 
significandy change between 10000 to 100000 simulation itera- 
tions (see Figure 8). The commitment of the phage to either lysis 
or lysogeny is assigned based on the relative abundance of cl and 
Cro proteins (> 1 for lysogeny, < 1 for lysis). 

The trajectory path of the phage is divided into a fine mesh. 
The direction of the next step is computed based on the 
combination of the deterministic force field and the stochastic 
noise. For every value of MOI, the ratio of the phages assigned to 
every outcome in the simulation is considered as the probability of 
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CI Concentration CI Concentration CI Concentration 

Figure 5. The effect of MOI on the Waddington landscape of phages in the lysis/lysogeny decision making. Three different values of 
MOI are presented in three rows. (A, B, C) The 2-D force-field representation for MOI =1,3 and 6 respectively. The concentration of cl and Cro free- 
monomers are assigned to the .Y-axis and j-axis respectively. The small arrows are directed toward the attractor states which are shown by the blue 
areas. (D, E, F) The 3-D potential function representation for MOI= 1, 3 and 6 respectively. The cl concentration, the value of the potential function, 
and the Cro concentration are shown by x-axis, ^-axis, and r-axis respectively. The attractor states are displayed as the local minimum of the 
potential function. 

doi:1 0.1 371 /journal.pone.01 03569.g005 



every decision (Figure 7). The results of the simulation show a 
strong agreement with both experimental results reported in [24] 
and our expectation from a rational player reported in Figure 4. 

To investigate the predictability of the model, the effect of the 
bacterium volume on the lysis/lysogeny decision is analyzed 
(Figure 7). In a population of bacteria with similar MOI, the 
computational model predicts that the smallest bacterium would 
have the highest probability of lysogeny. This is in agreement with 
in vitro experiments. It also shows that smallest volume bacterium, 
at the same MOI, has the highest viral concentration. It has a 
positive effect on cl expression and increases the chance of 
lysogeny [24]. 

Discussion 

We have studied the decision-making problem of bacteriophag- 
es lambda after infecting the E. coli bacterium. In particular, we 
have proposed a model of rational decision making that considers 



phages as competitive and cooperative agents and where there 
exists a trade-off between lysogenic and lytic pathways. The 
lysogenic pathway maintains the host bacterium as a host for 
reproduction, while lytic pathway increases the number of 
bacteriophages in the environment and destroys the host 
bacterium [22-24,28,29,37]. Thus, a rational decision maintains 
the host bacterium when there are enough bacteriophages in the 
environment and lysis otherwise. Making this decision requires 
external information from the environment which is represented 
by the average MOI. Such information is inferred based on 
internal signals such as size, MOI, and concentration of the host 
bacterium [38]. Given these parameters, we have analyzed the 
behavior of a rational decision in various environmental situations 
(see Figure 3). Our results show that the behavior of a rational 
player matches that of bacteriophages as reported in the 
experimental results of [24,30-32] (see Figure 4). 





Figure 6. The Waddington landscape of lysis/lysogeny in different MOI. The Waddington landscape changes from MOI= 1 to MOI = 7. 
The back-to-front axis represents increasing MOI values, the right-to-left axis shows different cl protein concentrations, and the vertical axis shows 
the pseudo-potential function. At MOI= 1 there is almost one valley with cl = 0 which represents a lysis attractor state. Increasing MOI results in a 
bifurcation and at MOI = 7 there are two distinct valleys with low (lysis) and high (lysogeny) cl concentrations. The potential values higher than —4 
are truncated for better visualization. 
doi:1 0.1 371 /journal.pone.01 03569.g006 
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Figure 7. The probability of lysogeny in different models. (A) Lysogeny decision-making probability based on Waddington's model affected 
by bacterium size. The horizontal axis is MOI value and the vertical axis shows the probability of lysogeny. The results are displayed for different £ 
coli sizes from 0.7, 0.9 and 1.1. (B, C) The probability of lysogeny when the size of the host bacterium is 0.9 and 1.1 respectively. Black curve 
represents the behavior of a rational player. Blue dashed curve shows the behavior of bacteriophages from the intracellular view based on its 
regulatory network. Red circles show the behavior of bacteriophages based on the experimental data of [24]. 
doi:1 0.1 371 /journal.pone.01 03569.g007 



From an intracellular point of view, the probability distribution 
of the lysis/lysogeny decision, under different MOI values, is 
defined by the structure of the gene regulatory network [30- 
35,39]. In this context, mapping the structure of the gene 
regulatory network to the Waddington landscape gives compre- 
hensive insight into the decision-making process. We trace how the 
environmental fluctuations reshape the landscape and change the 
fate of bacteriophages (see Figure 5 and Figure 6). We have 
argued that there exists an internal switch that makes decisions 
based on internal signals. We call it an internal biased coin that is 
encoded in the genome and inherited from its ancestors. 

Moreover, internal signals affect the weights of the edges in the 
gene regulatory network which, in turn, determine the shape of the 
Waddington landscape, the position of stable points, and the 
behavior of bacteriophages. On the other hand, internal signals 



also influence the estimation of average MOI which defines the 
expected utility function, and consequently, the mixed strategy of a 
rational player. Our results demonstrate that internal signals 
determine both the weights of edges in the regulatory network and 
the estimation of average MOI such that the probability 
distribution over attractors of the gene regulatory network is 
almost the same as the mixed strategy of a rational player (see 
Figure 7). 

Overall, the gene regulatory network controls the decision 
making of individuals in a biological environment [17,18]. Every 
cell makes decisions based on its intracellular genetic network and 
internal signals. The later determine the weights of the edges in 
gene regulatory network, and lead the organism to a stable point. 
A well-defined biochemical network, for decision making, is the 
one that can meet all of the cells requirements and changes the 
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Figure 8. Number of simulation runs for robust results. We perform measurement of lysis/lysogeny probabilities with 10' simulation runs 
where 0</<5. For each />0, we computed the average absolute change in measured probability values from 10 ,_1 to 10' simulation runs (x-axis 
labels show 10'). Our analysis shows 100000 simulation runs are sufficient for producing robust results. 
doi:1 0.1 371 /journal.pone.01 03569.g008 
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probability of their strategies based on different environmental 
situations. Mutations occur over time and change the structure of 
proteins, the shape of the regulatory network, the decision 
landscape, and the probability distribution of different decisions 
of an organism. Organisms with the highest adaptability will 
survive [19]. We have proposed a framework that shows how the 
gene regulatory network in the survivors leads the organism to 
behave almost the same as a rational player. 

Limitations of the Study, Open Questions, and 
Future Work 

The proposed framework can be used to study various 
fundamental cellular decision-making problems, e.g. decision 
between two different pathways of ATP synthesis in yeast. We 
believe that in binary decision-making problems, there is a bistable 
network which determines the fate of the organism. Internal 
signals define both the weights of the edges of the bistable network 
from intracellular perspective, and the estimation of the environ- 
mental state considering them as competitor and cooperative 
organisms. Based on the proposed framework, the bistable 
network can be viewed as a biased coin that defines an almost 
optimal mixed strategy in the presence of cooperators and 
competitors. More precisely, our approach can be used for other 
cellular decision-making processes to quantitatively assess how the 
corresponding gene regulatory network affects the rational 
decision-making process in different isogenic cells. 

We use sigmoid utility functions and the logit-response rule to 
model the rational player's decision-making process. Different 
utility functions and rules can be used for modeling decision- 
making processes. Studying the effect of different functions on 
decision-making models might be an appropriate direction for 
future studies. In addition, the gene regulatory dynamic equations 
presented in [22] are used as a basis for studying the decision- 
making process at the intracellular level, though other various 
mathematical models can also be used to represent the dynamic of 
a gene regulatory network in lysis/lysogeny decision making. 
Furthermore, the gene regulatory model we have used from [22] is 
limited to the triple key genes in lysis/lysogeny decision making. 
Investigating post-transcription modifications and translation 
noises that can cause deviations from quasi-steady state approx- 
imations of gene expression dynamics is a limitation of our study 
and can be studied further in the future. 

Conclusion 

In this paper we present a framework to show how the decision 
making of isogenic cells corresponds to the underlying gene 
regulatory network from an intracellular perspective. We use a 
game theoretic approach to model the decision-making process, 
considering cells as cooperative and competitive agents. We also 
quantitatively model the underlying gene regulatory network and 
use the Waddington landscape for a comprehensive understanding 
of attractor states. We demonstrate that the attractor states of the 
corresponding Waddington landscape fit the strategies of a 
rational player in the corresponding game. We study a 
fundamental decision-making problem in E. coli where a decision 
between lysis/lysogeny is made. We show that a decision between 
lysis or lysogeny from an intracellular perspective, considering the 
gene regulatory networks, almost matches a rational decision 
which maximizes the chance of living. To conclude, we provide a 
framework that uses game theory on one end and gene regulatory 
networks on the other end to study various fundamental cellular 
decision-making problems. The proposed framework explains how 



cells make decisions similar to a rational player in an analogous 
game. 

Materials and Methods 

Mathematical model of a rational decision 

To model the decision-making process of bacteriophages in E. 
coli, we first show how the environmental state can be estimated 
based on the internal information of phages such as MOI and size 
of the host bacterium. Then we present a method to determine the 
utility of each action and its corresponding probability based on 
the environmental state. 

Determining environmental state. A phage inside E. coli 
senses the size and the multiplicity of infection by probing the 
concentration of proteins and the concentration of dimers such as 
cl and Cro inside the host bacterium [33,34]. Therefore, in our 
model, size and MOI of the host bacterium are considered as 
internal signals. 

Let the random variable E be the average MOI within a 
population of phages and bacteria, where the average MOI is the 
total number of phages infecting a bacterium over the total 
number of bacteria in the environment. Also let II be the MOI of 
a host bacterium with size L, and [i be the number of phages 
inside the host bacteria. A bacterium with size L is a bacterium 
which its volume is L times the average volume of all bacteria in 
the environment A rational player infers about the average MOI 
and defines posterior probability f{E = X\Ii = p) based on the 
Bayes rule: 



f{E = X\I L = l i) = 



f(I L = n\E=X)f(E = X) 



(1) 



To compute the posterior probability, we estimate the value of 
f(lh = n\E ' = X),f(E ' = X), and /(II = Consider an environment 
with average MOI X in which phages move and infect bacteria 
randomly. The infection process can be seen as a Poisson process 
with average rate X [36] . The Poisson distribution is a probability 
distribution that expresses the probability of a given number of 
events that occur in a fixed interval of time and/ or space when the 
events occur with a known average rate [40] . On the other hand, 
the bacterium infection by phages is a random process, with an 
average number of infections A, which is independent of previous 
infections. Hence, the probability distribution of MOI for a single 
bacterium of size L, given the average MOI X, can be expressed by 
a Poisson distribution: 



f(I L = n\E = X) =/(/! =ji\E = LX) = 



\LXf 



(2) 



For prior probability f(E = X.) we assign a uniform distribution 
between 0 and A, i.e., f(E = X)= — for 0<X <A. The above 

simple distribution enjoys the maximum entropy and equivalendy 
the minimum knowledge among all priors. Uniform distribution is 
the most rational assumption when there is no prior information 
on the parameter [41,42]. At the end, probability /(II = can be 
obtained by marginalizing out X as follows: 



/(Il = H, E = x)dx - 



f(lL = fi\E = x)f (E = x)dx 
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(3) 



Now, probability f(E = X\lL = li) can be calculated based on 
equations 1, 2, and 3 as follows: 



f(E = l\I L =n)- 



f(I L = H\E = X)f(E = X) 

m=n) 



Le~ u (LXf 



(4) 



e x x tx dx 



From the definition of gamma function, we have 

rXu+l)= e~ x x fl dx, where r(fi+l) = [i\ for a positive integer 
Jo 

/•LA 

ji [43]. There is also a nice bound for the value of e~ x x^'dx as 

Jo 



follows: 



/•LA 

e- x x"dx=r( f i+ 1)-0((LA+ l)e- LA ). 

Jo 



(5) 



Equation 5 shows that T(^+l) is a good estimation for 

LA 

e~ x x >i dx. Putting all these facts together, we have: 



f(E=X\I L = n)*sLx 



e- L \L Xf 



LA' 



(6) 



(7) 



Utility Function. The lysis action increases the number of 
phages and kills the host bacterium. The lysogeny action only takes 
care of the host bacterium. Therefore, decision making in E. coli 
depends on the average MOI. It means if the number of phages in 
the environment is high and there are not enough host bacteria, a 
rational player takes care of host bacterium and lysogenize. 
Otherwise, a rational player tends to lyse. The utility function is 
modeled with a threshold a. The utility of lysis overcomes the 
utility of lysogeny if and only if a rational player infers that the 
average MOI is less than a.. 

To simplify the model, we look at the environmental state as a 
result of individual's strategies and define the utility functions 
based on the environmental state. We demonstrate the utility of 
individual j by function Uj : ,4x£— >R. Note that individual's 
strategies affect the state of the environment in the future. Thus, in 
our model the utility function implicitly depends on other 
individual's strategies. Let Uj(a, X) be the utility of action a e A 



of individual j £ Af when the environmental state, average MOI, is 
X. Let's assume the set of all possible actions is 
A = {lysis, lysogeny}. The utility of lysis should drop at threshold 
point a and the utility of lysogeny should rise at the same threshold 
point. This means the phages tend to keep the average MOI 
around a which is their desired average MOI. We focus on the 
sigmoid function and define the utility functions as follows: 



«,(lysis, X) = oP/(oP + X'), 



^■(lysogeny, X) = X" /(a" + 1"), 



(8) 



(9) 



where a is a threshold point and phages are more likely to choose 
lysogenic pathway when their estimation about the average MOI 
exceeds threshold a and rj is a parameter that defines the slope of 
the sigmoid function. 

Note that phages only detect internal signals, e.g., MOI and size 
of the host bacterium and thus they need an estimation about the 
utility of each action based on the observed internal signals rather 
than the average MOI. This is done by inferring about the average 
MOI based on the internal signals. Therefore, the expected utility 
of action a e {lysis, lysogeny} for a phage with host bacterium 
MOI \i and host bacterium size L, u*(a, Ii = p), is defined by 
integrating over all possible environmental average MOIs as 
follows: 



Uj(a, I L = n)=Ex[uj(a, X)\I L =n] = f(E = X\I L = n)Uj(a, X)dX. (10) 



Strategies. There are several rules in the population games 
and evolutionary dynamics that determine the strategy of players 
based on their utility. We employ the noisy best-response rule, the 
logit-response rule, as a well-known rule in the discrete choice 
literature for environmental evolution that well matches our 
setting [44-47]. In the logit-response rule, every individual plays 
its best-response strategy with a probability close to 1. However, 
we allow a small possibility for making mistakes. Individuals might 
make mistakes in their inferences if their information about their 
surroundings are noisy or the agents are not entirely rational 
[44,48-51]. In the logit-response rule the probability that 
individual j takes action a e {lysis, lysogeny}, Sj(a), is propor- 
tional to e H^^=^ as follows: 



Sj(a) -- 



Puj( a ,I L =n) 



^QymJ L =/i) 0KT(lysogeny,/ L =/i) ' 



(11) 



where /? e R + determines how noisy the system is. /?=oo shows 
that the system is noise-free and every individual plays its best 
action and /? = 0 represents a full noisy environment in which 
every individual plays randomly. A value between these two 
extreme points is chosen for modeling the behavior of real-world 
decision makers [44,48-51]. 

Mathematical model of a gene regulatory network 

The gene regulatory networks are mathematically represented 
in many different ways, including Boolean networks [52], Bayesian 
networks [53], ordinary differential equations (ODEs) [22], hybrid 
models [54], and even game theory [55]. Among all, we selected 
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ODEs since they represent the dynamical states of small networks 
more precisely along the time. 

The ordinary differential equations (ODEs) reported by [22] are 
used as a basis for the computational study of the lysis/lysogeny 
decision making at a molecular resolution. 

A non-restrictive quasi-steady state approximation (QSSA) 
method is applied to reduce the number of equations based on 
the great difference in the rate of fast (binding and unbinding of 
proteins) versus slow reactions (transcription, translation and 
protein degradation): 



dx 
dt 



dy 



fiax 



— bx, 



fiC 



dt \+x 2 +y : 



-dy. 



(12) 



Here n stands for multiplicity of infection (MOI), x and y are 
the rescaled concentrations of cl and Cro free monomers, a is 
transcriptional rate of cl when the promoter is bound, c is the 
transcriptional rate of Cro when the promoter is unbound, and b 
and d are degradation rates for cl and Cro proteins respectively. 
More details about equations and parameters are provided in the 
supplementing information. 

We model equation set 12 as a two-dimensional sample space, 
the force-field representation of which is provided in 5. Since this 
force-field representation is not a gradient vector fields, it is not 
possible to find an exact potential function for it. Therefore, we use 
a computationally feasible pseudo-potential function, only for a 
graphical representation of the landscape of the gene regulatory 
network. This potential function is monotonically increasing while 
it goes from an attractor center out to transient states (Figure 5). 
The logarithmic scale enables us to track the precise position of the 
attractor points. A pseudo-potential function is defined as follows: 



m.,)=io gU /(^ + (^ +( 



(13) 



Here e is a small positive value ensures that the logarithm 
function is computed for positive values. This equation assigns the 
minimum potential log(e) to all attractor states where dx/dt and 
dy I dt are zero. In theory, a similar potential function is assigned to 
unstable fixed points. However the vector-space representations of 
the force fields show such points are not present in our 
experiments (Figure 5, A-C). 

Additional Information 

Programs for modeling both a rational decision and gene 
regulatory network have been written in C++ and MATLAB. 



Supporting Information 

Figure SI The root-mean-square error. (A) This figure 
shows the root-mean-square error when /? = 3 and jj = 1.3 are 
fixed, and a varies from 0 to 5. (B) This figure shows the root- 
mean-square error when a = 2. 1 and )J = 1 .3 are fixed, and ft varies 
from 0 to 20. (C) This figure shows the root-mean-square error 
when a = 2.l and /? = 3 are fixed, and t] varies from 0.1 to 7. 
(TIFF) 

Figure S2 Lysogeny probability changes by different 
rate of cl dimerization. Decreased dimerization rate of cl will 
decrease the lysogeny probability. A factor of 0.3 to the cl 
dimerization rate decreases the depth of lysogenic attractor (left). 
By a factor of 0.1, the lysogenic valley vanishes and only the lytic 
attractor remains (right). 
(TIFF) 

Figure S3 Alternations to the network structure will 
destruct the function of the genetic switch. First, the 
negative effect of the Cro protein on the expression of its own gene 
is replaced by a positive effect. Second, the positive effect of cl on 
its own gene is replaced by a negative effect. 
(TIFF) 

Figure S4 Promoter mutations can alter the decision 
landscape. The mutations that increase the binding affinity of 
Cro to the operator sites will deepen the lytic attractor (left), while 
the contrary mutations increase the chance of lysogeny (right). 
(TIFF) 

Appendix SI Definition of parameters. In this appendix we 
define parameters for modeling a rational decision and the 
Waddington model. 
(PDF) 

Appendix S2 Verifying predictability of the Waddington 
model. In this appendix we perform several analyses for 
confirming the predictability of the proposed Waddington model. 
(PDF) 

Video SI The effect of MOI on the Waddington land- 
scape of phages in the lysis/lysogeny decision making. 

This video describes how the positions of attractor points change 
in the Waddinton landscape, for the host bacterium MOI from 1 
to 7. 
(MP4) 
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